/*******************************************************************************

	SSC0143 - Programacao Concorrente
	Prof. Julio Cesar Estrela
	
	T1 - Calculo de PI com 10000000 de casas decimais
	
	Grupo: 
	
	Daniel de Oliveira Pereira		6427171
	Guilherme Takasawa Yagui		6426698
	Vinicius Marques Stocco			6427233
	
	Março de 2013
	
********************************************************************************	
				Algoritmo de Borwein - Versao Concorrente	
*******************************************************************************/

#include <stdio.h>
#include <math.h>
#include <time.h>
#include <gmp.h>
#include <pthread.h>
#include <semaphore.h>

// numero de casas decimais corretas de PI
#define PD	10000000

// numero de bits minimo para armazenar PI com o numero de casas decimais corretas 
#define PB 	(int) ceil((log(10) / log(2)) * PD)

// numero de iteracoes necessario para calcular PI com o numero de casas decimais corretas desejado 
#define IT	(int) ceil(log(PD) / log(4))

// numero de iteracoes em espera
#define ITE 2

// tamanho maximo da fila
#define TMF 5

// posicao atual na fila
#define AT 	pos

// proxima posicao da fila
#define PR  (pos + 1) % TMF

mpf_t yk[TMF]; 					// fila para armazenar yk
mpf_t ak[TMF];					// fila para armazenar ak
mpf_t bk;						// pi

mpf_t tmp1, tmp2, tmp3, tmp4;	// auxiliares
	
pthread_mutex_t myk, mak;		// mutexes

sem_t sm, syk, sak;				// semafaros

pthread_t ty, ta;				// threads

unsigned int tfyk, tfak;		// tamanho atual da fila


// funcao para calcular y
// y(k+1) = (1-(1-y(k)^4)^(1/4))/(1+(1-y(k)^4)^(1/4))
void *Y(void *t)
{	
	unsigned int tfy2, tfy1 = 0;
	unsigned int pos = 0;
	
	while (1)
	{
		//printf("\n Y!");
		mpf_pow_ui(tmp1, yk[AT], 4);
		mpf_ui_sub(tmp1, 1, tmp1);
		mpf_sqrt(tmp1, tmp1);
		mpf_sqrt(tmp1, tmp1);
		mpf_set(tmp2, tmp1);
		mpf_ui_sub(tmp1, 1, tmp1);
		mpf_add_ui(tmp2, tmp2, 1);
		mpf_div(yk[PR], tmp1, tmp2);
		
		pos = PR;
		//tfy2 = tfy1;
		
		pthread_mutex_lock(&myk);
		tfy1 = ++tfyk;
		pthread_mutex_unlock(&myk);
		
		//printf("\t y = %d", tfy1);
		//if ((tfy1 == ITE) && ((tfy1 - tfy2) > 0)) 
		if (tfy1 == 1)	
			sem_post(&sak);
			
		if (tfy1 == TMF)
		{
			//printf("\t\t Yk cheio \n");
			sem_wait(&syk);
		}				
	}
}
  
// calcula a
// a(k+1) = ak*(1+y(k+1))^4-2^(2k+3)*y(k+1)*(1+y(k+1)+y(k+1)^2)     
void *A(void *t)
{
	unsigned int aux;
	unsigned int k = 0;
	unsigned int tfy;
	unsigned int tfa2, tfa1 = 0;
	unsigned int pos = 0;
	
	sem_wait(&sak);
	
	while (1)
	{
		//printf("\n A!");
		//aux 2 ^ (2k+3)
		aux = (unsigned int) pow(2.0, (2.0 * k + 3.0));
	
		mpf_add_ui(tmp3, yk[PR], 1);
		mpf_pow_ui(tmp3, tmp3, 4);
		mpf_mul(tmp3, tmp3, ak[AT]);
		mpf_pow_ui(tmp4, yk[PR], 2);
		mpf_add(tmp4, tmp4, yk[PR]);
		mpf_add_ui(tmp4, tmp4, 1);
		mpf_mul(tmp4, tmp4, yk[PR]);
		mpf_mul_ui(tmp4, tmp4, aux);
		mpf_sub(ak[PR], tmp3, tmp4);
		
		k++;
		pos = PR;
		//tfa2 = tfa1;
	
		pthread_mutex_lock(&mak);
		tfa1 = ++tfak;
		pthread_mutex_unlock(&mak);
		
		//printf("\t a = %d", tfa1);			
		//if ((tfa1 == ITE) && ((tfa1 - tfa2) > 0))
		if (tfa1 == 1)
			sem_post(&sm);
			
		pthread_mutex_lock(&myk);
		tfy = --tfyk;
		pthread_mutex_unlock(&myk);
		
		if (tfy == (TMF - 1))
			sem_post(&syk);
		
		//printf("\t yA = %d", tfy);
		
		if ((tfy == 0) || (tfa1 == TMF))
		{
			//printf("\t\t Yk vazio ou Ak cheio \n");
			sem_wait(&sak);		
		}	
	}
}

int main(char argc, char *argv[])
{
	FILE *f;
	
	clock_t t;		// tempo
	
	unsigned int i;
	unsigned int tfa;
	
	// define a precisão desejada para cada variavel
	mpf_set_default_prec(PB);
	
	// marca  o clock inicial da execucao do algoritmo
	t = clock();
	
	// cria o arquivo que ira armazenar os valores de PI
	f = fopen("borwein_paralelo.txt", "w");
	
	// aloca memoria para as variaveis utilizadas para calcular PI
	for (i = 0; i < TMF; i++)
	{
		mpf_init(ak[i]);
		mpf_init(yk[i]);
	}
	mpf_init(bk);
	mpf_init(tmp1);
	mpf_init(tmp2);
	mpf_init(tmp3);
	mpf_init(tmp4);
	
	// inicializa y(0)
	mpf_sqrt_ui(tmp1, 2);
	mpf_sub_ui(yk[0], tmp1, 1);
	
	// inicializa a(0)
	mpf_sqrt_ui(tmp2, 2);
	mpf_mul_ui(tmp2, tmp2, 4);
	mpf_ui_sub(ak[0], 6, tmp2);
	
	// tamanho atual das filas de yk e ak
	tfyk = 0;
	tfak = 0;
	
	//iniciliaza os mutexes
	pthread_mutex_init(&myk, NULL);
	pthread_mutex_init(&mak, NULL);
	
	// inicializa os semaforos
	sem_init(&sm, 0, 0);
	sem_init(&syk, 0, 0);
	sem_init(&sak, 0, 0);
	
	// cria as threads utilizadas para calcular y(k+1) e a(k+1)
	pthread_create(&ta, NULL, &A, NULL);
	pthread_create(&ty, NULL, &Y, NULL);
	
	// bloqueia a execução da main
	sem_wait(&sm);
	
	for (i = 0; i < IT; i++)
	{
		//printf("\n M!");
		// calcula PI = b(k) = 1/a(k+1)
		mpf_ui_div(bk, 1, ak[(i+1) % TMF]);
		
		pthread_mutex_lock(&mak);
		tfa = --tfak;
		pthread_mutex_unlock(&mak);
		
		//printf("\t aM = %d", tfa);
		
		if (tfa == (TMF - 1))
			sem_post(&sak);
		
		if (tfa == 0)
		{ 
			//printf("\t\t Ak vazio \n");
			sem_wait(&sm);	
		}
		//gmp_printf("\n %.*Ff", PD, bk);
		gmp_fprintf(f, "%.*Ff \n", PD, bk);
	}
	printf("\n Numero de iteracoes : %d \n", IT);
	
	// cancela a execucao das threads
	pthread_cancel(ty);
	pthread_cancel(ta);
	
	// libera recursos alocados pelas threads
	pthread_join(ty, NULL);
	pthread_join(ta, NULL);
	
	// destroi os semaforos
	sem_destroy(&sm);
	sem_destroy(&syk);
	sem_destroy(&sak);
	
	//destroi os mutexes
	pthread_mutex_destroy(&myk);
	pthread_mutex_destroy(&mak);
	
	// libera a memoria alocada para calcular PI
	for (i = 0; i < TMF; i++)
	{
		mpf_clear(ak[i]);
		mpf_clear(yk[i]);
	}
	mpf_clear(bk);
	mpf_clear(tmp1);
	mpf_clear(tmp2);
	mpf_clear(tmp3);
	mpf_clear(tmp4);
	
	// fecha o arquivo com os valores de PI
	fclose(f);
	
	// calcula o numero clocks necessario para executar o algoritmo
	t = clock() - t; 
	
	printf("\n Tempo de execucao : %.3f s \n\n", (double) t / CLOCKS_PER_SEC);

	return 0;
}
